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Abstract 

Wc consider a size-structured model for cell division and address the question of determining 
the division (birth) rate from the measured stable size distribution of the population. We propose a 
new regularization technique based on a filtering approach. We prove convergence of the algorithm 
and validate the theoretical results by implementing numerical simulations, based on classical tech- 
niques. We compare the results for direct and inverse problems, for the filtering method and for 
the quasi-reversibility method proposed in [1]. 

1 Introduction 

The use of size-structured models to describe biological systems has attracted the interest of many 
authors and has a long standing tradition. In particular, the use of size structures was very well 
documented and compared to experiments in the 70's. This led to the survey book [2] and subsequent 
mathematical analysis (see also the references in [3]). Needless to say, in such models it is crucial for 
the analysis, computer simulation and prediction to calibrate the corresponding model parameters so 
as to obtain good quantitative results. Indeed, in the inverse problem literature, a number of authors 
have addressed the calibration of certain structured population models. See for example [H O O [7] 
and references therein. 

In this article, we consider theoretical and numerical aspects of the inverse problem of determining 
the division rate coefficient B = B[x) in the following specific size-structured model for cell division: 

^n(t, x) + -i^n{t, x) + B{x)n{t, x) = 4B{2x)n{t, 2x), x ^ 0, t ^ 0, 
n{t,x = 0) = 0,t>0, (1) 
n{0,x) = n^{x) > 0. 
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Here, the cell density is represented by n{t, x) at time t and size x. The division rate B expresses the 
division of cells of size 2x into two cells of size x. 

By making use of flux cytometry technologies for instance, it is possible to determine cell populations 
with certain properties as protein content on a large scale of tenths of thousands of cells. In other 
applications, like coagulation fragmentation equation [HI El [lOl [HI [12], or prion aggregation and 
fragmentation |13lll4t [T5]. similar equations arise, and much less is known on aggregate size repartition. 
The division rate B{x), on the contrary, is not directly measurable. 

The long time behavior of solutions is well known. Indeed, it was proved in |161I17| that under fairly 
general conditions on the coefficients, there is a unique solution (A^, Aq) to the following eigenvalue 
problem 



^N + {Xo + B{x))N = AB{2x)N{2x), x ^ 0, 



N{x = 0) = 0, 

N{x) > for X > 0, 

where Aq > and iVe^"" e L°° n L'^ for all fJ. < Xq. 
It was shown in [18\ 116] 

n{t,x)e~^°^ moN{x), in L^(R+,(j){x)dx), 

where the weight (j) is the unique solution to the adjoint problem 



(2) 



^'^ + (Ao + 5(x))0 = 25(x)</.(^ 



dx 

(x) > 0, 



(t){x)N{x)dx 



2^" 

1. 



X ^ 0, 



(3) 



In other words, Aq is the growth rate of such a system and is usually called "Malthus parameter" in 
population biology. From [181 [El [3] we also know that Aq is related to by the relation 



An 



r Ndx _ 

Jq°° xNdx 



(4) 



The question we address here is the following: How can we estimate the division rate B from the 
knowledge of the steady dynamics N and Aq ? The inverse problem thus consists of finding B a 
solution to 

d 

ABi2x)N(2x) - B(x)N(x) = Lix) := ^N(x) + AoiV(x), x ^ 0, (5) 

ox 

assuming that (A^, Aq) is known, or thanks to ((H) that is known. As seen in [1], this problem is 
well-posed if A^ satisfies strong regularity properties such as -^^ix) £ L^(M+) for some p ^ 1. 

However, in practical applications we have only an approximate knowledge of {N, Aq), given by noisy 
data (A'g, A^), with G L^(R+) for instance. This means that we have no way of controlling -^N^, 
so we cannot control the precision of a solution i?^ to problem ([5]) when a perturbed Nf, replaces A^. 
Furthermore, it is not even clear whether such a B^ exists. 



^Actually, our knowledge of Ao is presumably an order of precision higher than that of A*', since the rate Ao can be 
estimated independently by means of time information. 
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The question we focus on is then: How to approximate the problem ([5]) in order to get a solution 
as close as possible to the exact division rate B? 



We remark that, in the context of noisy data, the inverse problem under consideration is ill-posed [T] 
and thus regularization would be required. A natural tool to be invoked from the inverse problem 
literature would be some kind of Tikhonov regularization method |191l20j. However, this would lead to 
computationally intensive problems. Indeed, for each forward problem evaluation a dilation-differential 
equation of the form ([2]) would have to be solved. 

In [l], two of the present authors proposed a method of regularization consisting in the solution of 
the following approximate problem: 

a^(i?.,aiVe) +4B.,a(y)iV.(2/) = i?.,a(f)iV.(f) + AoA^ed) +2^(^Ar,(|)), y > 0, 

{Be,aNem = 0, 

where a is a regularizing parameter. It was shown that a convergence rate of order ^/e could be 
obtained, for a = 0{y/e), where e is the error on the data in an appropriate norm. 

The above method of the solution to the inverse problem will be called quasi reversibility in accor- 
dance with the general spirit of the terminology of [21^ [22] . The main goal of this work is to investigate 
the numerics of such approach, to consider an alternative technique based on filtering ideas and to 
compare the performance of the different methods. The alternative technique is also analyzed from 
the theoretical point of view and estimates are presented. 

In this work, we have modified slightly the original regularization equation by writing Ag^Q, instead 
of Aq for the reasons we shall explain in the sequel. Thus, we work with 



a^(Se,aiV,) + iB,MNe{y) = A,(|) + Ae,„A,(|) + 2 9- (^iV,(|)) , y > 0, 



(6) 



^ {Be,aN,){0) = 0. 



Indeed, in order to conserve regularity properties of the solution H = BN to the inverse problem, we 
want it to be both in L^(M+) and in L^(Kj^,xdx) in order to express that both the total number of 
cells and the total biomass are finite. Hence, formal integration of Equation ([6]) gives 

Xe,a / Nedx = / Be,aNedx, (7) 

Jo Jo 
and integration against the weight x gives 

poo poo poo 

-a Be^aNedx = 4Ae,„ / xNedx -4 Nedx. (8) 
Jo ' Jo Jo 

Hence, we have to choose, according to the eigenvalue theory: 

f°° N dr 

''^'^ xN^dx + N^dx ^""^ 

The choice of A^^q can be understood as a compatibility condition when a > and for q = it tells us 
that {N, Aq) is overdetermined data for the inverse problem. Therefore, if we have a priori knowledge 
on Aq, we could verify its distance to A^ q, as a way of checking the error of the inverse problem solution. 
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The plan of this work is the following: In Section [2l we propose yet another method to regularize 
the inverse problem, and obtain a convergence rate. The convergence rate turns out to be as good as 
the one in [1]. In Section [3] we give a numerical method to solve it, and in Section U] we show some 
numerical simulations so as to compare the accuracy of the different methods. 



2 Regularization by Filtering 

2.1 Filtering approach 

Taking a closer look at Equation ([5]), we see that all the difficulties come from the differential term 
-^N. In [U, the choice was to add an equivalent derivative a-^i^B to the equation; here on the 
contrary, we choose to regularize it by a convolution method. 
For a > 0, we use the notation 

p^{x) = -p{-), peC^im, / p{x)dx = l, p^O, Supp(/>) C [0,1], (10) 
and we replace in ([5]) the term ^-/Vg + XqNi; by 

{g^^e + Xe,aNe) * Pa{x) = * {g^Pa + K,aPa){x) = J Ns{x'){—pa + Xe,aPa){x - x')dx' . 

We now use the notation 

Ne,o, = Ne* Pa- 
in this way, we obtain a smooth term in L^(]R+). Furthermore, N^^a converges to Ne in L^(IR+) when 
a tends to zero. We now have to consider the following problem: 
Find -Bg^Q solution of 

ABe,ai2x)Ne,ai2x) + (x)iV,,a (x) = ^iV,,„ + A,,„7V,,„ (x) , X ^ . (11) 

As in Equation ([6]), for the quasi-reversibility method, we need to choose Xs^a appropriately. Indeed, 
we perform the same manipulations leading to Equation ([9]) to get 

!^ X^s,a{x)dx 

By Theorem IA.3I (see the Appendix), we know that the problem in Equation ()lip has a unique 
solution Be,a e L'^{'^+,Nl^dx). 

2.2 Estimates for the filtering approach 

The main result of this section establishes an estimate for the regularization of the inverse problem 
by means of the filtering method described above. 

Theorem 2.1 Suppose that N G H'^{^+) and B e L~(M+), B ^ verify (0). Let e > and 
Ne G -L^(M+), Ne{x) >Oforx>0, such that 

\\Ne - N\\l2^^^) ^ e\\N\\L2 
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Let Be,a G L'^(R+,N^adx) be the unique solution of and We have the following estimate: 

C 

\\Be,a - B\\L2(^Nldx) ^ C{a + |Ae,Q - Ao| ) | |iV| |h2(R+) + ~ ^IIl2(R+)i (13) 

where C is a constant depending only on \\B\\loo^ \\Be,a\\L°° and the regularizing function p. 
This theorem rehes on a first estimate. 

Proposition 2.2 Using the same notations as in Theorem we have 

m^aN,,a-BN\\l,^,^^^C{l + Xl) + \\Ne-N\\l,^,^)+C{a^+\>^e,a-Xo\')\\N\\l,^^^) , (14) 
where C depends only on the regularizing function p. 

Proof of Prop. [23} Denote hy Q = Bs,aNe,a - BN, R = N^^a - N and 5 = X^^a - Xq. From 
Equations ^ and pT|) . Q verifies: 

£r{x) + XoR{x) + SN,,a{x) + Q{x) = 4Q{2x), x ^ 0, 

(15) 

Q{x = 0) = 0. 

(Since N^^a £ H^{R^), the definition of Q{x = 0) is not ambiguous.) Multiplying this equation by 
Q{2x) and integrating on the interval (0, y) yields 

y y y 

4 J Q{2xfdx = j Q{2x)-^R{x)dx + Xo J Q{2x)R{x)dx 



y 



+ 6 J Q{2x)Ne,a{x)dx + j Q{2x)Q{x)dx. 





From the Cauchy-Schwarz inequality, after the change of variables j; — > 2x, we have 
y y r, y y y 



Ao f Qi2x) 



4 J Q{2xfdx (^^) (^)'^^ + \j Qi'^xfdx + yJ CR{xfdx + ^ J ^ 

GO 

y 

„ y y y 2 



^2 



2 



j Ne,a{xfdx + ]^j Q{2xfdx + ^j Q{2xfdx + j Q{2xfdx. 



00 
We take, for instance, C = Xq. We obtain 

\\Be,aNe,a " BN\\l, ^ ||iV, * —p^ - q^N\\1, + Xl\\Ne * Pa " N\\l, + |A,,, - Xo\''\\Ne * PaWh ■ (16) 

The last two terms of this inequahty are easy to estimate, writing 

\\N,*Pa-N\\L2 ^ \\Ne*Pa-N*Pa\\L2 + \\N*Pa-N\\L2^C{\\N,-N\\L2+a\\N\\Hi), 



and 

\\N,*p^\\l2^C\\N\\l2. 
It remains to evaluate the first term on the right-hand side of inequality (|16p . We write 



By a convolution estimate we evaluate the first term as 

Since \-§^Pa{x)\dx = ^ \-^p{y)\dy, we have 

d d C(p) d d 

To evaluate the last term ||A^* -g^Pa — ^-^11^2; we extend to M the functions N and by zero and 
consider their Fourier transforms. We denote /(^) the Fourier transform of / G L^(M+) at i^, where / 
is extended as zero on M_. We obtain by Fourier analysis 

Using that 

l^^^^^2^KC(p)a^ (17) 
where C{p) only depends on the regularization function p, we have that 

Going back to p^ . this concludes the proof of Proposition 12.21 □ 
We can now deduce the proof of Theorem 12.11 We write: 

\\Be,a - ^||L2(^r|rf^.) ^ \\Be^a]^e - Ns^a\\L'2{R+) + 1 1 -Be,a A''e,a - B N\\i^2q^^^ + \\B N - B iVg 1 1 ^2 . 

Using Proposition 12.21 and the fact that 

\\Ne - iVe,a||L2 ^ 2||iV, - iV| + a| | 1 1 j^i , 

this inequality gives the result. □ 

3 Numerical Solution of the Inverse Problem 

This section is concerned with the numerical aspects of the solution of the inverse problem. In order 
to do that we start with a description of the solution to the direct one in Subsection 13.11 
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3.1 Direct Problem 



In the direct problem, we assume we know the prohferation rate B, we look for N and Aq > solutions 
of ([2|). For this purpose, we solve the time-dependent problem ([T|) and look for a steady dynamics. As 
already said, this problem is well-posed (see for instance [3]) and it was proved in [18] that solutions 
grow at an exponential rate towards pN{x)e^"^ with p = n{0,x)(l){x)dx, recalling the notation in 
([3]). Furthermore, under more restrictive conditions it was shown in p!6] that there exists constants 
/X > and C(n°) > 0, such that 

||n(t,a;)e-^°* -/>iV(x)||ii(R+,^(,)d,) ^ Ce"^'*. 

To solve it numerically, we discretize the problem ([1]) along a regular grid, denote by At the time step 
and by Ax = L/I the spatial step, where / denotes the number of points and L the computational 
domain length: Xi = iAx, ^ i ^ I. 



We use an upwind finite volume method (cf. [231 [2 

^»+^ At 
'^i = j '^ik'^t^y)dy, j n{kAt + s,x-^i)ds 



X 1 



For the time discretization, we use a marching technique. We choose the time step At so as to satisfy 
the largest possible CFL stability criteria '■= = ^■ 
The numerical scheme is given, for i = 1, /, by ng = and 



i + ' + = i?2i-in^,„i + 2B2,ni + i32^+lr^^.+l , (18) 

At Ax 

with the convention that rij = for j > /. For stability reasons, we have used an implicit method for 
the division term in the left hand side and explicit for the right hand side of the equation. The specific 
form for the right hand side is simply motivated by the need of also dividing cells of odd labels. 
According to the power algorithm, we do not keep n^^^ from (llSp but rather renormalize it as 



~k+i 



n 



I 

Ax ^ 



It is standard, for these positive matrices arising in (jlSp . that 

/ 

n 



^=+1 ^ iV, 2^iVi = l, A^^>0, 



i=l 



where is the dominant eigenvector for the problem 

A^^ - Ni.i 



+ (Ao + Bi)N, = B2^-lN2i-l + 2B2iN2^ + B2^+lN2i+l. 



Ax 

One can also find the dominant eigenvalue as 



i=l 
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For matrices with one dominant eigenvalue and a corresponding one-dimensional eigenspace, it is 
known that the power algorithm is fast and in fact converges with exponential rate |26j . In practice 
we can stop the iterations when the relative error on the normalized quantity 




is small enough, say of the order of 10 



3.2 Inverse Problem: General Strategy 

In the sequel, we denote by H the product B.N and its approximations. Indeed, from Equations ([6]) 
or ([TT]) . we have to search for the product H = Bg^^aNe or H = Bf^^a^e^a before computing B^^^- In 
particular, we cannot avoid a loss of information where Nf, is small, i.e., for x ~ or x S> 1. 
The inverse problem ([5]), as well as (jlip . can be written as 

4H{2x) - H{x) = L{x), (19) 

with different expressions for H and L. We may think of two possible numerical approaches. 

Strategy 1. Compute H{2x) from H(x): This means that we re-write Equation ()19p with the new 
variable y = 2x, and arrive at 

iH{y)-H{^) = L{^). (20) 

The scheme departs from zero, and one deduces the values of Hi step by step, from the knowledge of 
Hj for j ^ i — 1. 

Strategy 2. Compute H{x) from H{2x): The scheme departs from the largest point x = L of our 
simulation domain. We suppose that for x ^ L we have H{x) = H{L) = (it is relevant since we 
suppose that N vanishes for x large: see below), and then deduce the smaller values Hi step by step, 
from the knowledge of Hj for j ^ i + 1. 



The two approaches do not necessarily lead to the same result because the continuous equation 

4H{2x) - H{x) = (21) 

has infinitely many solutions. This issue is interesting on its own and is related to the construction of 
wavelets, see |27j . It is discussed in Proposition I A . 1 1 of the Appendix. 

By imposing H E L^(M+), we select a unique solution, as shown in Theorem IA.3[ The question is 
then: Which numerical strategy should we use to select the correct solution, i.e. the one in L^(R-|_) ? 

Among the solutions of Equation (jlOp . we single out two, defined by the power series: 

+ 00 +00 

i/a)(2;) = ^2"2"L(2""x) and = _^22"L(2"x) , V x > 0. 

71=1 71=0 

Proposition lA.ll shows that for L G L^(R+, x^dx), there is a unique solution in L^(M_|_, x^dx), given 
by if p < 3 and by i?^^) if 

p > 3 (and the power series converge in the corresponding spaces). 
For i? > smooth and bounded from above and from below, we know that N is smooth and 
vanishes at x ~ and x ~ oo, and BN inherits these properties. For instance, we know that 
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H G Lp'{dx) n L^[x^dx). By uniqueness of a solution in each space, Proposition lA.ll implies that 
= H^^\ or equivalently: 

+00 

J2 22"L(2"x) = 0, Vx^O. 

n=—oo 

This very particular property cannot be verified at the discrete level. Hence, the two strategies 
generally give two different approximations of the same solution of (jl9p . The first strategy selects an 
approximation of the solution H^^^ whereas the second selects an approximation of the solution H^"^^ . 
In the case of a very regular data N, then if(^) will perform better around infinity, whereas H^^^ will 
be better around zero. However, if is a solution of Equation when we increase the number of 
points, the two approaches converge to the same solution since 

Since our simulation domain [0, L] is bounded and contains zero, we prefer the first strategy. This 
choice is confirmed by all the numerical tests we have performed: the second approach has always 
lead to a solution exploding around zero. However, for the sake of completeness, we also describe the 
scheme we used for the second approach. 



3.3 Inverse Problem: Filtering Approach 

According to strategies 1 and 2, we now present two approaches to handle the numerical solution of the 
inverse problem regularized with the filtering approach. Both need to first compute the convolution 
terms arising in (jlip . To do so we first take the Fast Fourier Transform F of Ag, multiply it by i^PaiO, 
and then take the inverse Fast Fourier Transform F* . We choose and define the regularization function 
Pa by its Fourier transform: 

This leads us to the numerical approximation 

d_ 

dx 

We also impose dN^fi = for compatibility with the continuous equation and further use. 

As mentioned earlier, there are two alternatives, either starting from zero or coming from infinity. 



Ne,a ~ dNa = F* ( i^PaimNem ] ■ (22) 



The Filtering Approach Starting from Zero (strategy 1). We solve Equation (jlip considered 
as an equation in the variable y = 2x, that is to say ()20p . in order to compute its solution H^^\x). 
At the discrete level, we use the notations 



The discrete version of ([20]) reads 

AhI = h{+l{, VO^i^/, (23) 



and we need to define the quantities Gi . We choose 

2 



G± when i is even, 

G, = { ' , (24) 

2 (Gi-i + Gi+i j when z is odd. 



2 
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In particular, we have Hq = 0. 

Summing up all the terms in ()23p for 1 ^ i ^ /, we find (with / even to simplify): 



1=0 i=0 



Since we have assumed that has exponential decay for x ^ 1, it follows that 

^h( = ^l{ + Ej, with \Ei\^2Y^\Hi\. (25) 



i=0 i=l ,_i 



Multiplying (j23p by and summing up again, we find 

f / 
Y,XiL{ = Fi, with IF/K (26) 

i=l ,-I 

As a consequence, we can choose: 

as the discrete version of the relations or (jl2p . 

The Filtering Approach Starting from Infinity (strategy 2). Another method is to discretize 
the formulation (jl9p in order to compute its solution H^'^\x). We define the extension = for 
i ^ I + 1, and for 2 < i < /, we define by backward iterations 

f f f f 

This however does not apply to the indices i = 0, 1 and we set = -j- = and H( = AHr^ — L'[ . 



By summing up all the terms in (I28p . we find balance properties equivalent to (I25p - (l26p . but with 
remainders Ei and Fj depending on Hi and H2 instead of H^i . One has to check a posteriori that 
these last quantities are very small ; it is not the case in a standard calculation, but becomes true 
when the precision of the direct problem scheme increases. 

3.4 Inverse problem: Quasi-Reversibility Approach 

In this section, we present a numerical scheme for the regularized inverse problem proposed in [T]. 
This problem leads to solving ([6]) taken at y = 2x, that is 

a^(S,,,iV,) + iB,MNe{y) = B.,4f )iV,(f ) + A,,,iV,(f ) + 2^ (iV,(f )) , y > 0, 

where a > is the regularizing parameter and A^^q is defined by Q. This gives, in a discretized 
version, after dropping the index e, 

= ^^^^ftjENl- ^^^^ 
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For the numerical discretization we set H^^ = and also recall that A'^o = ^-nd assume that the 
data satisfies Nj+i = 0. We use a standard upwind scheme for the differential term: 

« (^Q - H^_^) + 4H^ = Hf + L?, (30) 

Ax 2 2 



where we have defined the fractional indices as in the filtering approach by ()24p . and here 



Ax 



If we neglect the terms , we can easily verify a discrete version of the balance laws and ([9]), 

equivalent to (f25]) -(l26l). 



4 Numerical Tests 

As input data, we take the values of the function N obtained by the numerical solution of the direct 
problem in Section 13. H we add a random noise uniformly distributed in [—§,§], and we enforce 
nonnegativity of the data 

Ne = max(iV + er, 0). 

We solve the direct problem on a regular grid of / + 1 points, on an interval [0, 2L]. We need L large 
enough, such that it is possible to assume that A^(x ^ L) ~ and we have checked it a posteriori. 
Indeed, we have seen that this property is essential when we use the inverse schemes on a domain 
[0, L] in order to verify the balance laws ([Tj)-®. In other words, we solve the direct problem on a 
domain twice larger than for the inverse problem. In the numerical tests we take L = 4, and we show 
the numerical solution N only on the interval [0, L] since it is uniformly small on [L, 2L]. 

We solve the inverse problem by the different methods on a regular grid of /i + 1 points on [0, L], 
with Axi = L/Ii. This grid is taken ten times finer than the grid used for the direct problem, i.e. we 
take Ii = 101. Since we have chosen L large enough so that N(x ^ L) ~ 0, we have always obtained 
that indeed H{x ^ L) ^ 0. 

As before, we denote by and the solution data H obtained respectively by the quasi- 
reversibility method of Section 13.41 and by the first filtering approach (from zero) of Section 13. 3[ We 
also define a solution H-^^ by mixing both methods, i.e. by solving the following equation: 



* Pa{x), X ^ 0, 

(31) 



B,^^{x = 0)N,{x = 0) = 0, 
where Xs^a is defined by 



^^^^ C * pgdx ^ ^^^^ 

xNs * Padx + f /g~ Ne * Padx ' 



The relative error is measured, as seen in Theorem 12.11 and in Theorem 5.1 of [T], by 

Q ^ \\BN,-HQ\\f, ^ \\BN,-Hf\\l ^ \\BN,-Hf^ 
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We have divided by ||A'e||j;^2 and not by ||A^||j|^2 because in practice we only know the entry data with 
noise. 

In order to illustrate the accuracy of our method, we also compare it to a naive way (brute force) of 
considering the equation. Namely, we approximate -^N{x) by a second-order Euler scheme without 
regularization. It gives a solution by the same formula (j30p . where we simply take a = 0. 



The Direct Problem. We have first tested the direct problem for various division rates B. Three 
different solutions for three given division rates B are depicted in Figure [T] with 800 grid points. 

In the particular case when B is constant, we can go further and evaluate the computational error. 
Then, we know that X = B and the exact solution A/'oxact can be explicitly calculated, as shown in 
[m E] ) by the formula: 

oo 

iVexact(x)=ivJ^a„e-2"^^ (33) 

n=0 

where the coefficients are defined recursively by = 1 ^md a„ = (— l)"^Srx"' ^ chosen to 
ensure the mass one normalization. We take B = 1 and obtain the continuous curve of Figure [TJ We 
can measure here the relative error by 

ll-Aoxactll/i 

where N represents the numerical solution of Section [XTl We choose this norm because for B constant, 
the solution of the adjoint problem is = 1 and the General Relative Entropy Principle (|181 13]) gives 
us that this quantity decreases along the time iterations. Still for 800 points, we obtain 6^ = 7.7.10"^. 




NforB=1 

N for B=1 for )(<15, then linear increase to B=5 
NforB=1+exp(-8(x-2A 



5 1 Is 2 2.5 3 3.5 4 




Figure 1: Solutions (left) obtained by the numerical resolution of Section [3. II for the direct problem 
with three different division rates B (right). 



The noiseless case (e = 0). In the simplest case where the data is perfectly known, i.e. for e = 0, 
we verify that the different schemes allow us to recover B. Since the precision of the data is directly 
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linked to the number of points used in the scheme, we run the codes with 1.000 points for the direct 
problem (below, we will take only 100 points). 

We test several values of a and we use the three functions B of Figure [U for each method for the 
inverse problem. The error estimate is found to depend on the method used but not significantly on 
the division rate B. Therefore we have drawn in Figure [2] the average error estimates for the three 
division rates B. In Figure [3] we have depicted the products B.N in the case B = 1 and a = 0.01 
(other cases are similar): it shows that the precision obtained is satisfactory. In Figures HJ [5] and [6] we 
have drawn the approximations of B in each of the three cases, calculated only for N > 0.01 (indeed, 
for too small the division leads to insignificant results on B). 

Not surprisingly, the brute force method reveals to be satisfactory, with an error estimate of 6^ = 
1.3.10"^, since we are in the case where N is very regular. The filtering method can reach this level of 
error for a = 10~^ but cannot go further. However, both the quasi reversibility method and the mixed 
method given by Equation (j3ip improve it with minimum values 6^ = 6.9.10"^ and 6-1^^ = 6.5.10"^ 
reached for a = 10~^. 



-6- Quasi-reversibility method 

Fitering method 
-ir- Filtering then quasi-reversibility method 

Bruit force 




-0.2 



■— Quasi-Reversibility method. a=0 01, error 5=3.10'^ 

— Filtering approach. o:=0 05. Error 5=9.10'^ 

— Filter (o.=0.05)-t Quasi-Reversibility (a=0.01). Error 6=8,10"- 
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Figure 2: For e = 0, numerical errors obtained 
with the different methods for the inverse prob- 
lem. 



Figure 3: Numerical reconstruction of B.N ob- 
tained by each method for the inverse problem 
when B = 1, s = and a = 0.01. 



Link between the noise level e and the regularization parameter a. For noise levels e = 0.01, 
e = 0.05 and e = 0.1 respectively, the Figures [71 [8] and [9] give the curves e as a function of a for the 
three inverse methods. We compare the reconstructed division rates B in Figures [TUl and [TTl 

Each of the error curves presents a minimum for an optimal value of a, as expressed by estimate 
()13p for instance. In Figures [T2l [T3l [H] and \T5\ we have compared three curves, drawn in a log-log 
scale: ^/e to serve as a reference curve, /(e) = min5(a,e), and g{e) = argmin (5(a, e). One can see 

a a 

that for each method, these three curves have comparable slopes (^ on a log-log scale): they show that 
even though the combination of filtering and quasi-reversibility method improves the optimal errors 
in absolute value, it does not change the order of convergence of the approximation, which remains 
of order 0{^/e). Figure [TCI gives also the convergence of the filtering method for much smaller values 
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Exact problem: B=1 

B^^, Quasi-Reversibility method 

B^^, Filtering approach 

B ,Fiiter + Quasi-Reversibility 




■ Exact problem. B=1 then 5 
^, Quasi-Reversibility method 

. ^. Filtering approach 

Filter + Quasi-Reversibility 



DO 




Figure 4: Reconstructed division rate B using 
the three inverse methods, for e = 0, a = 0.01 
with computed from B = I. 



Figure 5: Reconstructed division rate B, for 
e = 0, a = 0.01 and a jump S = 1 to 5 as in 
Figure [T] 



of e (for which an increased number of 500 points has been taken, in order to avoid numerical bias): 
the comparison with ^/e is there particularly evident, and we have obtained similar curves for the 
two other methods. The speed of convergence is though in complete accordance with the theoretical 
estimate (fT3l) . 

Influence of the choice of Aq instead of Xe^a- To evaluate the influence of the error term due to 
the distance {X^^a — Xo\, we compare the curves obtained respectively by taking on the inverse code 
the exact Aq or the value Xg^a expressed by the balance laws. They are drawn in Figure [13] for the 
quasi-reversibility method. They show that even though the a priori knowledge of Aq improves the 
error in absolute value, it does not change the order of convergence of the scheme. Thus it is in 
complete accordance with Estimate (fT3]) . 

5 Conclusions 

We have considered size-structured equations connected to several areas of biology from cell division 
to prion proliferation by aggregation and fragmentation. We have addressed the numerical efficiency 
of some inverse problem solution methods to tackle the problem of recovering the division rate from 
the size distribution of cells. The latter involves a dilation equation with a singular right-hand side 
that needs regularization for actual implementation. For that purpose, we have introduced a filtering 
method and proved its convergence for noisy data. This method brings in an operator that has a 
non-trivial kernel and we have selected a numerical approximation that is able to recover the natural 
solution we want to reach. 

The implementation of the inverse algorithm, based on the filtering method, confirms the conver- 
gence analysis. In particular, there is an optimal regularization parameter as can be seen in the 
graphs of Figures [71 [8] or [9] for instance. Comparison with a quasi-reversibility method introduced 
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evolution of error with a, e=0.01 
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6=0 . 


■ — Exact problem; B=exp(-8[x-2)'^)+l 
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Figure 6: Reconstructed division rate B, for 
£ = 0,a = 0.01 and B = l + exp{-8{x - 2)2). 



Figure 7: Numerical error when e = 0.01 for 
the different methods. 



earher leads to the conclusion that a combination of filtering and quasi-reversibility methods seems to 
be more efficient because the oscillations are reduced, but without improving the rate of convergence. 

We also analyzed the impact of using the exact value of Aq or the on the different solutions of the 
inverse problem. In our simulations, the difference between using Aq or Ag seemed to be immaterial as 
far as the accuracy of method is concerned. This is in perfect accordance with the theoretical estimate 

m- 

The above remarks open several directions for continuation and extension of the present work. On 
the practical side, the present work sets the stage for the use of experimental data either from the 
existing literature or from more recent biological experiments. On the theoretical side, the possibility 
of improving the convergence by combining the filtering and quasi-reversibility methods should be 
investigated further. 

Finally, we point out that although the Tikhonov method is more standard, we did not study it 
so far because it seems more time consuming. Indeed, iterations are needed to solve both the direct 
problem and the inverse one. To overcome such difficulty a completely new theory has to be developed 
so as to suit the particular structure of our model. This provides yet another direction for future work. 
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Figure 8: Numerical error when e = 0.05. Figure 9: Numerical error when e = 0.1. 

A Well-Posedness of Functional Equation Associated to the Inverse 
Problem. 

We have seen that the regularization method for the inverse problem relies mostly on solving the 
equation 

AH{2x) - H{x) = L{x), X > 0. (34) 

Eventhough this equation is formally very simple, its analysis reveals some complexity. It may admit 
several solutions in general. Among them, we can mention two with simple representation formulas 
(we leave to the reader to check they are indeed formally solutions) 

(x) = ^ 2-2"L(2-"x). (35) 

n=l 

+ 00 

/7(2) (^) = _ ^ 22"L(2"x), (36) 

n=0 

To clarify this issue and motivate our choice of a solution, we first state general results concerning 
solutions to (j34p and then come back to our original problem ([5]). 
We first mention the following 

Proposition A.l Let L € (IR+ , x^dx) , with p ^ 3, then there exists a unique solution H G 
L^{^+,xPdx) to ^ and 

• for p < 3, this solution is given explicitly by the formula Ii35\) . Furthermore, for I < q < oo, if 
L e Li{R+) then H^^^ G Li{R+). 

• for p > 3, this solution is given explicitly by the formula ^36\). 

Because we look for an integrable function H (the number of cells is supposedly finite), the function 
H^^^ is preferable (take q = 1). It also behaves better near rr ~ because the weight p < 3 imposes 
that H^^^ vanishes at as we expect. 
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E=0.01,B=exp(.8(x-2)^)+1 




Figure 10: In the case e = 0.01, a = 0.05, B = 1 + e numerical solution B.N (left) and B 

(right) by the different methods. 



From the point of view of exact solutions of the direct problem, we find that H = BN and L belong 
to all spaces L^(]R_|., x^dx) for all p G M. Therefore, the two solutions coincide and in principle we 
could choose any of them. In practice, errors on the data L are better handled by //^^^ than by //^^^ 
for the afore mentioned reason. Notice indeed that these two solutions are different in general. One 
can check for instance that for L = 0, there is a singular distributional solution 5^=0 • Furthermore, 

Lemma A.2 The solutions to with L = in P'(0, oo) have the form with f £ P'(M) a 

log(2)— periodic distribution. 

Proof of Proposition IA.lt We consider the Hilbert space X = L^(R^_, x^dx) and we simply apply 
the Lax-Milgram theorem to a properly chosen bilinear form. 



Case 1, p < 3. We solve the equation in the variable y = 2x that is (j20p . and consider the bi- 
linear form a(u, v) on X x X defined by 

rco roc y 

a{u,v) = A u{y)v{y)yPdy - u{-)v{y)y^dy. 







'2' 

This form is obviously continuous and it remains to prove that it is coercive. We have 



oo 



a(n,n)=4/ uiyfyUy - i tx(^)n(y)rt ^ (4 - 2^) / uiyfy^dy, 
Jo Jo ^ Jo 

p+1 

and it is indeed coercive as long as a = 4 — 2 2 is positive which holds true for p < 3. The Lax- 
Milgram Theorem asserts that there is a unique H £ X such that a{H, •) = (L, •), where (•, ••) denotes 
the inner product in X, that is a solution of ()20p . 

Case 2, p > 3. We work in the variable x and consider the continuous bilinear form 6(n, v) on 
X X X defined by 

poo poo 

b{u,v) = —4 / u{2x)v{x)x^dx + / u{x)v{x)x^dx. 
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the case e 





and B = 1, the numerical solu- 
tion B.N by the different methods. 

The same calculation leads us to: 

j-OO 

b{u,u) u{xf'x^dx, 
Jo 

and the same conclusion holds. 



Filtering + quasi-reversibility method 




Figure 12: Filtering and quasi-reversibility 
method: minimal error as a function of the 
noise level e for the optimal value a, with a 
comparison to the theoretical curve ^/e. 



with a = 1 - 2 2 > 0, 



To check formulae (j36p and (j35p . it remains to prove that these solutions belong to the corresponding 
spaces: 

oo oo 
n=l n=l 

This sum converges iff p > 3. In the same way, we write: 

oo oo 
n=0 n=0 

which converges ijf p > 3. □ 

Proof of Lemma IA.2t When L = 0, we first define 7i G T^'{0, oo) as the second antiderivative of H, 
and notice that it should verify 

n{2x) = n{x). 

We perform the change of variables y = log(2;) and notice that, if W G P'(0, oo), it is equivalent to 
look for solutions / E D'(M) of 

/(y + log(2)) =/(y). (37) 
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Figure 13: Quasi-reversibility method, with (left) A^^q given by relation ([9]) or (right) a priori knowl- 
edge of Ao : minimal error and optimal regularization parameter a as functions of the noise level e, 
with a comparison to the theoretical curve ^/e. We see that the a priori knowledge of Aq does not 
improve the speed of convergence of the scheme. 



Hence, all the solutions in X''(0, oo) are given by —^—^ — -, where / G P'(]R). □ 

To conclude this Appendix, we come back to our original problem ^ and draw the consequences 
in terms of B, not H. 



Theorem A. 3 Let N £ L^{R+), with N{x) > for x > 0. Let L £ L^{R+). There exists a unique 
B e L'^(R+,N'^dx) solution of 



Proof: The theorem follows directly from Proposition lA.il for p = 0, and since > 0, we can define 
B = H/N for B G L'^{R+, N^dx). □ 

This theorem shows that we can find a solution B of ^ for all N and all A, this is the basis of our 
algorithm. However, if we want that the solution B belongs to the space -L"^(R+; xN{x)dx), integration 
of ([38]) multiplied by x shows that L has to satisfy the condition 



Applying this to Equation we recover that Aq = /q°° N{x)dx/ xN{x)dx. In the case of Equa- 
tions dl]) and (jlip respectively, we get formulae ([9]) and (I12|) . which discrete versions are expressed by 
([29]) and 

In view of these considerations, it is better to use a discrete scheme defined by a matrix A that 
preserves a similar discrete property. Namely, for all H = (Hi), we should have '^i{AH)i = 0, in 



other words the vector of components i belongs to the kernel of the adjoint of A. Indeed, this property 
yields the (discrete) regularity H G L^ {M.^; xdx) . 




AB{2x)N{2x) - B{x)N{x) 



L{x). 



(38) 
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Figure 14: Filtering method for standard levels Figure 15: Filtering method for smaller values 

of noise. of e and increased number of points. 
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